MODULE STATISTICS
  ! SALVADOR NAVARRO-LOZANO
  ! JANUARY 26, 2004
  IMPLICIT NONE
!!$  PRIVATE :: P_Sort_Stat,P_Inssort_Stat ! Helper function, same as in matrix module
  
  INTERFACE Mean
     MODULE PROCEDURE MEAN_V, MEAN_M
  END INTERFACE
  
  INTERFACE Variance
     MODULE PROCEDURE VARIANCE_V, VARIANCE_M
  END INTERFACE
  
  INTERFACE Standard_Deviation
     MODULE PROCEDURE STDEV_V, STDEV_M
  END INTERFACE
  
CONTAINS
  
  REAL(8) FUNCTION MEAN_V(V)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: V(:)
    MEAN_V = SUM(V) / DBLE(SIZE(V))
  END FUNCTION MEAN_V
  
  FUNCTION MEAN_M(M)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: M(:,:)
    REAL(8) :: MEAN_M(SIZE(M,2)),n
    INTEGER :: j
    n=DBLE(SIZE(M,1))
    DO j = 1, SIZE(M,2)
       MEAN_M(j) = SUM(M(:,j)) / n
    END DO
  END FUNCTION MEAN_M
  
  REAL(8) FUNCTION VARIANCE_V(V)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: V(:)
    REAL(8) :: m
    m = MEAN(V)
    VARIANCE_V = (SUM((V-m)*(V-m))) / (DBLE(SIZE(V))-1.0d0)
  END FUNCTION VARIANCE_V
  
  FUNCTION VARIANCE_M(M)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: M(:,:)
    REAL(8) :: VARIANCE_M(SIZE(M,2))
    INTEGER :: j
    DO j = 1, SIZE(M,2)
       VARIANCE_M(j) = VARIANCE_V(M(:,j))
    END DO
  END FUNCTION VARIANCE_M
  
  REAL(8) FUNCTION STDEV_V(V)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: V(:)
    STDEV_V = SQRT(VARIANCE_V(V))
  END FUNCTION STDEV_V
  
  FUNCTION STDEV_M(M)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: M(:,:)
    REAL(8) :: STDEV_M(SIZE(M,2))
    STDEV_M = SQRT(VARIANCE_M(M))
  END FUNCTION STDEV_M

  REAL(8) FUNCTION Covariance(V1,V2)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: V1(:),V2(:)
    REAL(8) :: m1,m2,a1(SIZE(V1))
    m1 = MEAN(V1)
    m2 = MEAN(V2)
    a1 = (V1 - m1)*(V2-m2)
    Covariance = SUM(a1) / (SIZE(V1)-1.0d0)
  END FUNCTION Covariance
  
  REAL(8) FUNCTION Correlation(V1,V2)
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: V1(:),V2(:)
    REAL(8) :: m1,m2,a1
    a1 = Covariance(V1,V2)
    m1 = Standard_Deviation(V1)
    m2 = Standard_Deviation(V2)
    Correlation = a1/(m1*m2)
  END FUNCTION Correlation

  FUNCTION Variance_Covariance(A) RESULT(B)
    IMPLICIT NONE 
    REAL(8), INTENT(IN) :: A(:,:)
    REAL(8) :: B(SIZE(A,2),SIZE(A,2))
    INTEGER :: i,j
    DO i = 1, SIZE(A,2)
       DO j = 1, SIZE(A,2) 
          B(i,j) = Covariance(a(:,i),a(:,j))
       END DO
    END DO
  END FUNCTION Variance_Covariance

  FUNCTION Histogram(vector,limit)
    ! H=HIST(vector,limit) returns the histogram of vector
    ! where H(1) returns the number of elements of vector <= limit(1)
    ! where H(i) returns the number of elements of vector between limit(i-1) and limit(i) (closed) for i<=SIZE(limitin)
    ! and H(SIZE(limit)+1) returns the number of elements of vector>limit(SIZE(limit))
    ! !!!!!!!!!!!Notice: H is an integer!!!!!!!!!!!!!!!!!!!!!!!
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: vector(:),limit(:)
    INTEGER :: Histogram(SIZE(limit)+1)
    INTEGER :: i,n
    n=SIZE(limit)
    Histogram=0
    Histogram(1)=COUNT(vector<=limit(1))
    DO i = 2, n
       Histogram(i)=COUNT((vector>limit(i-1)).AND.(vector<=limit(i)))
    END DO
    Histogram(n+1)=COUNT(vector>limit(n))
  END FUNCTION Histogram

  FUNCTION Centile(x,p)
    ! calculate the percentiles you ask for in p of data vector x
    ! p then is an real(8) vector between 1 and 100 (so you can ask for the 12.5 centile if you want)
    ! centile, pin must be of the same size. Uses the same method as Matlab 6.1
    IMPLICIT NONE
    REAL(8), INTENT(IN) :: x(:),p(:)
    REAL(8) :: centile(SIZE(p)),pp(SIZE(p))
    REAL(8) :: xx(SIZE(x)+2),q(SIZE(x)+2),step,h(SIZE(x)+1),dxx(SIZE(x)+1)
    REAL(8) :: s(SIZE(p)),del(SIZE(x)+1)
    INTEGER :: nx,nxx,np,i,ind(MAX(SIZE(x),SIZE(p))),k(SIZE(p))
    np=SIZE(p)
    nx=SIZE(x)
    nxx=nx+2
    step = 100.0d0/DBLE(nx)
    q(1) = 0.0d0
    q(nxx) = 100.0d0
    q(2) = 0.5d0*step
    DO i = 3, nx+1
       q(i) = q(i-1) + step
    END DO
    CALL P_Sort_Stat(x,xx(2:nx+1),ind(1:nx))
    xx(1)=xx(2)
    xx(nxx)=xx(nxx-1)
    DO i = 1, nx+1
       h(i) = q(i+1) - q(i)
       dxx(i) = xx(i+1) - xx(i)
    END DO
    ind=(/(i,i=1,SIZE(ind))/)
    CALL P_Sort_Stat(p,pp,ind(1:np))
    k = 0
    DO i = 2, nxx
       WHERE ((pp>=q(i-1)).AND.(pp<q(i))) k=i-1
    END DO
    WHERE (pp==q(nxx)) k=nxx-1
    s = pp - q(k)
    del = dxx/h
    centile(ind(1:np)) = xx(k) + s*del(k)
  END FUNCTION Centile

!!$  function dinvnr ( p, q )
!!$
!!$!*****************************************************************************80
!!$!
!!$!! DINVNR computes the inverse of the normal distribution.
!!$!
!!$!  Discussion:
!!$!
!!$!    This routine returns X such that 
!!$!
!!$!      CUMNOR(X) = P, 
!!$!
!!$!    that is, so that 
!!$!
!!$!      P = integral ( -Infinity <= T <= X ) exp(-U*U/2)/sqrt(2*PI) dU
!!$!
!!$!    The rational function on page 95 of Kennedy and Gentle is used as a 
!!$!    starting value for the Newton method of finding roots.
!!$!
!!$!  Reference:
!!$!
!!$!    William Kennedy, James Gentle, 
!!$!    Statistical Computing,
!!$!    Marcel Dekker, NY, 1980,
!!$!    QA276.4 K46
!!$!
!!$!  Parameters:
!!$!
!!$!    Input, real ( kind = 8 ) P, Q, the probability, and the complementary
!!$!    probability.
!!$!
!!$!    Output, real ( kind = 8 ) DINVNR, the argument X for which the
!!$!    Normal CDF has the value P.
!!$!
!!$  implicit none
!!$
!!$  real(8) :: ccum
!!$  real(8) :: cum
!!$  real(8) :: dinvnr
!!$  real(8) :: dx
!!$  real(8), parameter :: eps = 1.0D-13
!!$  integer :: i
!!$  integer, parameter :: maxit = 100
!!$  real(8) :: p
!!$  real(8) :: pp
!!$  real(8) :: q
!!$  real(8), parameter :: r2pi = 0.3989422804014326D+00
!!$  real(8) :: strtx
!!$  real(8) :: stvaln
!!$  real(8) :: xcur
!!$
!!$  pp = min ( p, q )
!!$  strtx = stvaln ( pp )
!!$  xcur = strtx
!!$!
!!$!  Newton iterations.
!!$!
!!$  do i = 1, maxit
!!$
!!$    call cumnor ( xcur, cum, ccum )
!!$    dx = ( cum - pp ) / ( r2pi * exp ( -0.5D+00 * xcur * xcur ) )
!!$    xcur = xcur - dx
!!$
!!$    if ( abs ( dx / xcur ) < eps ) then
!!$      if ( p <= q ) then
!!$        dinvnr = xcur
!!$      else
!!$        dinvnr = -xcur
!!$      end if
!!$      return
!!$    end if
!!$
!!$  end do
!!$
!!$  if ( p <= q ) then
!!$    dinvnr = strtx
!!$  else
!!$    dinvnr = -strtx
!!$  end if
!!$
!!$  return
!!$end function dinvnr
!!$
!!$function eval_pol ( a, n, x )
!!$
!!$!*****************************************************************************80
!!$!
!!$!! EVAL_POL evaluates a polynomial at X.
!!$!
!!$!  Discussion:
!!$!
!!$!    EVAL_POL = A(0) + A(1)*X + ... + A(N)*X**N
!!$!
!!$!  Modified:
!!$!
!!$!    15 December 1999
!!$!
!!$!  Parameters:
!!$!
!!$!    Input, real ( kind = 8 ) A(0:N), coefficients of the polynomial.
!!$!
!!$!    Input, integer N, length of A.
!!$!
!!$!    Input, real ( kind = 8 ) X, the point at which the polynomial 
!!$!    is to be evaluated.
!!$!
!!$!    Output, real ( kind = 8 ) EVAL_POL, the value of the polynomial at X.
!!$!
!!$  implicit none
!!$
!!$  integer n
!!$
!!$  real(8) :: a(0:n)
!!$  real(8) :: eval_pol
!!$  integer :: i
!!$  real(8) :: term
!!$  real(8) :: x
!!$
!!$  term = a(n)
!!$  do i = n - 1, 0, -1
!!$    term = term * x + a(i)
!!$  end do
!!$
!!$  eval_pol = term
!!$
!!$  return
!!$end function eval_pol
!!$
!!$function stvaln ( p )
!!$
!!$!*****************************************************************************80
!!$!
!!$!! STVALN provides starting values for the inverse of the normal distribution.
!!$!
!!$!  Discussion:
!!$!
!!$!    The routine returns an X for which it is approximately true that 
!!$!      P = CUMNOR(X),  
!!$!    that is, 
!!$!      P = Integral ( -infinity < U <= X ) exp(-U*U/2)/sqrt(2*PI) dU.
!!$!
!!$!  Reference:
!!$!
!!$!    William Kennedy, James Gentle,
!!$!    Statistical Computing, 
!!$!    Marcel Dekker, NY, 1980, page 95,
!!$!    QA276.4 K46
!!$!
!!$!  Parameters:
!!$!
!!$!    Input, real ( kind = 8 ) P, the probability whose normal deviate 
!!$!    is sought.
!!$!
!!$!    Output, real ( kind = 8 ) STVALN, the normal deviate whose probability
!!$!    is approximately P.
!!$!
!!$  implicit none
!!$
!!$  real(8) :: eval_pol
!!$  real(8) :: p
!!$  real(8) :: sgn
!!$  real(8) :: stvaln
!!$  real(8), parameter, dimension(0:4) :: xden = (/ &
!!$    0.993484626060D-01, &
!!$    0.588581570495D+00, &
!!$    0.531103462366D+00, &
!!$    0.103537752850D+00, &
!!$    0.38560700634D-02 /)
!!$  real(8), parameter, dimension(0:4) :: xnum = (/ &
!!$    -0.322232431088D+00, &
!!$    -1.000000000000D+00, &
!!$    -0.342242088547D+00, &
!!$    -0.204231210245D-01, &
!!$    -0.453642210148D-04 /)
!!$  real(8) :: y
!!$  real(8) :: z
!!$
!!$  if ( p <= 0.5D+00 ) then
!!$
!!$    sgn = -1.0D+00
!!$    z = p
!!$
!!$  else
!!$
!!$    sgn = 1.0D+00
!!$    z = 1.0D+00 - p
!!$
!!$  end if
!!$
!!$  y = sqrt ( -2.0D+00 * log ( z ) )
!!$  stvaln = y + eval_pol ( xnum, 4, y ) / eval_pol ( xden, 4, y )
!!$  stvaln = sgn * stvaln
!!$
!!$  return
!!$end function stvaln

!!!!!!!!!!!!!!!!!!!!!! HELPER FUNCTIONS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE P_Sort_Stat(Xin,xout,ind)
    !  Sorts xout into ascENDing order - _quicksort_stat
    REAL (8), DIMENSION (:), INTENT (in) :: Xin
    INTEGER, INTENT(INOUT) :: ind(:)
    REAL (8), DIMENSION (SIZE(ind)), INTENT (Out) :: xout
    xout=XIN
    Call P_SUBSOR_Stat(xout,1,Size(xout),ind)
    Call P_Inssort_Stat(xout,ind)
  CONTAINS
    RECURSIVE SUBROUTINE P_SUBSOR_Stat (xout, IDEB1, IFIN1,ind)
      !  Sorts xout from IDEB1 to IFIN1
      REAL(8), DIMENSION (:), INTENT (InOut) :: xout
      INTEGER, INTENT (In) :: IDEB1, IFIN1
      INTEGER, INTENT(INOUT) :: ind(:)
      INTEGER, PARAMETER :: NINS = 16 ! Max for insertion sort
      INTEGER :: ICRS, IDEB, IDCR, IFIN, IMIL,IWRK,IPIV
      REAL(8) :: XPIV, XWRK
      IDEB = IDEB1
      IFIN = IFIN1
      !  IF we DOn't have enough values to make it worth while, we leave
      !  them unsorted, and the final insertion sort will take care of them
      IF ((IFIN - IDEB) > NINS) THEN
         IMIL = (IDEB+IFIN) / 2
         !  One chooses a pivot, median of 1st, last, and middle values
         IF (xout(IMIL) < xout(IDEB)) THEN
            XWRK = xout (IDEB)
            IWRK = IND (IDEB)
            xout (IDEB) = xout (IMIL)
            IND (IDEB) = IND (IMIL)
            xout (IMIL) = XWRK
            IND (IMIL) = IWRK
         END IF
         IF (xout(IMIL) > xout(IFIN)) Then
            XWRK = xout (IFIN)
            IWRK = IND(IFIN)
            xout (IFIN) = xout (IMIL)
            IND (IFIN) = IND (IMIL)
            xout (IMIL) = XWRK
            IND (IMIL) = IWRK
            IF (xout(IMIL) < xout(IDEB)) Then
               XWRK = xout (IDEB)
               IWRK = IND (IDEB)
               xout (IDEB) = xout (IMIL)
               IND (IDEB) = IND (IMIL)
               xout (IMIL) = XWRK
               IND (IMIL) = IWRK
            END IF
         END IF
         XPIV = xout (IMIL)
         IPIV = IND (IMIL)
         !  One exchanges values to put those > pivot in the END and
         !  those <= pivot at the beginning
         ICRS = IDEB
         IDCR = IFIN
         ECH2: DO
            DO
               ICRS = ICRS + 1
               IF (ICRS >= IDCR) Then
                  !  the first  >  pivot is IDCR
                  !  the last   <= pivot is ICRS-1
                  !  Note: IF one arrives here on the first iteration, then
                  !  the pivot is the maximum of the set, the last value is equal
                  !  to it, and one can reduce by one the size of the set to process,
                  !  as IF xout (IFIN) > XPIV
                  Exit ECH2
               END IF
               IF (xout(ICRS) > XPIV) Exit
            END DO
            DO
               IF (xout(IDCR) <= XPIV) Exit
               IDCR = IDCR - 1
               IF (ICRS >= IDCR) Then
                  !  The last value < pivot is always ICRS-1
                  Exit ECH2
               END IF
            END DO
            XWRK = xout (IDCR)
            IWRK = IND (IDCR)
            xout (IDCR) = xout (ICRS)
            IND (IDCR) = IND (ICRS)
            xout (ICRS) = XWRK
            IND (ICRS) = IWRK
         END DO ECH2
         !  One now sorts each of the two sub-intervals
         Call P_SUBSOR_Stat (xout, IDEB1, ICRS-1,IND)
         Call P_SUBSOR_Stat (xout, IDCR, IFIN1,IND)
      END IF
      RETURN
    END SUBROUTINE P_SUBSOR_STAT
  END SUBROUTINE P_Sort_Stat

  SUBROUTINE P_Inssort_Stat (xout,IND)
    !  Sorts xout into increasing order (Insertion sort)
    REAL(8), DIMENSION (:), INTENT (InOut) :: xout
    INTEGER, INTENT(INOUT) :: ind(:)
    INTEGER :: ICRS, IDCR,IWRK
    REAL(8) :: XWRK
    DO ICRS = 2, Size (xout)
       XWRK = xout (ICRS)
       IWRK = IND (ICRS)
       IF (XWRK >= xout(ICRS-1)) Cycle
       xout (ICRS) = xout (ICRS-1)
       IND (ICRS) = IND (ICRS-1)
       DO IDCR = ICRS - 2, 1, - 1
          IF (XWRK >= xout(IDCR)) Exit
          xout (IDCR+1) = xout (IDCR)
          IND (IDCR+1) = IND (IDCR)
       END DO
       xout (IDCR+1) = XWRK
       IND (IDCR+1) = IWRK
    END DO
    RETURN
  END SUBROUTINE P_Inssort_Stat

  SUBROUTINE median(x, n, xmed)

    ! Find the median of X(1), ... , X(N), using as much of the quicksort
    ! algorithm as is needed to isolate it.
    ! N.B. On exit, the array X is partially ordered.
    
    !     Latest revision - 26 November 1996
    IMPLICIT NONE
    
    INTEGER, INTENT(IN)                :: n
    REAL, INTENT(IN OUT), DIMENSION(:) :: x
    REAL, INTENT(OUT)                  :: xmed
    
    ! Local variables
    
    REAL    :: temp, xhi, xlo, xmax, xmin
    LOGICAL :: odd
    INTEGER :: hi, lo, nby2, nby2p1, mid, i, j, k
    
    nby2 = n / 2
    nby2p1 = nby2 + 1
    odd = .true.
    
    !     HI & LO are position limits encompassing the median.
    
    IF (n == 2 * nby2) odd = .false.
    lo = 1
    hi = n
    IF (n < 3) THEN
       IF (n < 1) THEN
          xmed = 0.0
          RETURN
       END IF
       xmed = x(1)
       IF (n == 1) RETURN
       xmed = 0.5*(xmed + x(2))
       RETURN
    END IF
    
    !     Find median of 1st, middle & last values.
    
10  mid = (lo + hi)/2
    xmed = x(mid)
    xlo = x(lo)
    xhi = x(hi)
    IF (xhi < xlo) THEN          ! Swap xhi & xlo
       temp = xhi
       xhi = xlo
       xlo = temp
    END IF
    IF (xmed > xhi) THEN
       xmed = xhi
    ELSE IF (xmed < xlo) THEN
       xmed = xlo
    END IF
    
    ! The basic quicksort algorithm to move all values <= the sort key (XMED)
    ! to the left-hand end, and all higher values to the other end.
    
    i = lo
    j = hi
50  DO
       IF (x(i) >= xmed) EXIT
       i = i + 1
    END DO
    DO
       IF (x(j) <= xmed) EXIT
       j = j - 1
    END DO
    IF (i < j) THEN
       temp = x(i)
       x(i) = x(j)
       x(j) = temp
       i = i + 1
       j = j - 1
       
       !     Decide which half the median is in.
       
       IF (i <= j) GO TO 50
    END IF
    
    IF (.NOT. odd) THEN
       IF (j == nby2 .AND. i == nby2p1) GO TO 130
       IF (j < nby2) lo = i
       IF (i > nby2p1) hi = j
       IF (i /= j) GO TO 100
       IF (i == nby2) lo = nby2
       IF (j == nby2p1) hi = nby2p1
    ELSE
       IF (j < nby2p1) lo = i
       IF (i > nby2p1) hi = j
       IF (i /= j) GO TO 100
       
       ! Test whether median has been isolated.
       
       IF (i == nby2p1) RETURN
    END IF
100 IF (lo < hi - 1) GO TO 10
    
    IF (.NOT. odd) THEN
       xmed = 0.5*(x(nby2) + x(nby2p1))
       RETURN
    END IF
    temp = x(lo)
    IF (temp > x(hi)) THEN
       x(lo) = x(hi)
       x(hi) = temp
    END IF
    xmed = x(nby2p1)
    RETURN
    
    ! Special case, N even, J = N/2 & I = J + 1, so the median is
    ! between the two halves of the series.   Find max. of the first
    ! half & min. of the second half, then average.
    
130 xmax = x(1)
    DO k = lo, j
       xmax = MAX(xmax, x(k))
    END DO
    xmin = x(n)
    DO k = i, hi
       xmin = MIN(xmin, x(k))
    END DO
    xmed = 0.5*(xmin + xmax)
    
    RETURN
  END SUBROUTINE median

subroutine cumnor ( arg, cum, ccum )

!*****************************************************************************80
!
!! CUMNOR computes the cumulative normal distribution.
!
!  Discussion:
!
!    This function evaluates the normal distribution function:
!
!                              / x
!                     1       |       -t*t/2
!          P(x) = ----------- |      e       dt
!                 sqrt(2 pi)  |
!                             /-oo
!
!    This transportable program uses rational functions that
!    theoretically approximate the normal distribution function to
!    at least 18 significant decimal digits.  The accuracy achieved
!    depends on the arithmetic system, the compiler, the intrinsic
!    functions, and proper selection of the machine dependent
!    constants.
!
!  Author: 
!
!    William Cody
!    Mathematics and Computer Science Division
!    Argonne National Laboratory
!    Argonne, IL 60439
!
!  Reference:
!
!    William Cody,
!    Rational Chebyshev approximations for the error function,
!    Mathematics of Computation, 
!    1969, pages 631-637.
!
!    William Cody, 
!    Algorithm 715: 
!    SPECFUN - A Portable FORTRAN Package of Special Function Routines 
!    and Test Drivers,
!    ACM Transactions on Mathematical Software,
!    Volume 19, Number 1, 1993, pages 22-32.
!
!  Parameters:
!
!    Input, real(8) ARG, the upper limit of integration.
!
!    Output, real(8) CUM, CCUM, the Normal density CDF and
!    complementary CDF.
!
!  Local Parameters:
!
!    Local, real ( kind = 8 ) EPS, the argument below which anorm(x) 
!    may be represented by 0.5 and above which  x*x  will not underflow.
!    A conservative value is the largest machine number X
!    such that   1.0D+00 + X = 1.0D+00   to machine precision.
!
  implicit none

  real(8), parameter, dimension ( 5 ) :: a = (/ &
    2.2352520354606839287D+00, &
    1.6102823106855587881D+02, &
    1.0676894854603709582D+03, &
    1.8154981253343561249D+04, &
    6.5682337918207449113D-02 /)
  real(8) :: arg
  real(8), parameter, dimension ( 4 ) :: b = (/ &
    4.7202581904688241870D+01, &
    9.7609855173777669322D+02, &
    1.0260932208618978205D+04, &
    4.5507789335026729956D+04 /)
  real(8), parameter, dimension ( 9 ) :: c = (/ &
    3.9894151208813466764D-01, &
    8.8831497943883759412D+00, &
    9.3506656132177855979D+01, &
    5.9727027639480026226D+02, &
    2.4945375852903726711D+03, &
    6.8481904505362823326D+03, &
    1.1602651437647350124D+04, &
    9.8427148383839780218D+03, &
    1.0765576773720192317D-08 /)
  real(8) :: ccum
  real(8) :: cum
  real(8), parameter, dimension ( 8 ) :: d = (/ &
    2.2266688044328115691D+01, &
    2.3538790178262499861D+02, &
    1.5193775994075548050D+03, &
    6.4855582982667607550D+03, &
    1.8615571640885098091D+04, &
    3.4900952721145977266D+04, &
    3.8912003286093271411D+04, &
    1.9685429676859990727D+04 /)
  real(8) :: del
  real(8) :: eps
  integer :: i
  real(8), parameter, dimension ( 6 ) :: p = (/ &
    2.1589853405795699D-01, &
    1.274011611602473639D-01, &
    2.2235277870649807D-02, &
    1.421619193227893466D-03, &
    2.9112874951168792D-05, &
    2.307344176494017303D-02 /)
  real(8), parameter, dimension ( 5 ) :: q = (/ &
    1.28426009614491121D+00, &
    4.68238212480865118D-01, &
    6.59881378689285515D-02, &
    3.78239633202758244D-03, &
    7.29751555083966205D-05 /)
  real(8), parameter :: root32 = 5.656854248D+00
  real(8), parameter :: sixten = 16.0D+00
  real(8) :: temp
  real(8), parameter :: sqrpi = 3.9894228040143267794D-01
  real(8), parameter :: thrsh = 0.66291D+00
  real(8) :: x
  real(8) :: xden
  real(8) :: xnum
  real(8) :: y
  real(8) :: xsq
!
!  Machine dependent constants
!
  eps = epsilon ( 1.0D+00 ) * 0.5D+00

  x = arg
  y = abs ( x )

  if ( y <= thrsh ) then
!
!  Evaluate  anorm  for  |X| <= 0.66291
!
    if ( eps < y ) then
      xsq = x * x
    else
      xsq = 0.0D+00
    end if

    xnum = a(5) * xsq
    xden = xsq
    do i = 1, 3
      xnum = ( xnum + a(i) ) * xsq
      xden = ( xden + b(i) ) * xsq
    end do
    cum = x * ( xnum + a(4) ) / ( xden + b(4) )
    temp = cum
    cum = 0.5D+00 + temp
    ccum = 0.5D+00 - temp
!
!  Evaluate ANORM for 0.66291 <= |X| <= sqrt(32)
!
  else if ( y <= root32 ) then

    xnum = c(9) * y
    xden = y
    do i = 1, 7
      xnum = ( xnum + c(i) ) * y
      xden = ( xden + d(i) ) * y
    end do
    cum = ( xnum + c(8) ) / ( xden + d(8) )
    xsq = aint ( y * sixten ) / sixten
    del = ( y - xsq ) * ( y + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) * exp ( -del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if
!
!  Evaluate ANORM for sqrt(32) < |X|.
!
  else

    cum = 0.0D+00
    xsq = 1.0D+00 / ( x * x )
    xnum = p(6) * xsq
    xden = xsq
    do i = 1, 4
      xnum = ( xnum + p(i) ) * xsq
      xden = ( xden + q(i) ) * xsq
    end do

    cum = xsq * ( xnum + p(5) ) / ( xden + q(5) )
    cum = ( sqrpi - cum ) / y
    xsq = aint ( x * sixten ) / sixten
    del = ( x - xsq ) * ( x + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) &
      * exp ( - del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if

  end if

  if ( cum < tiny ( cum ) ) then
    cum = 0.0D+00
  end if

  if ( ccum < tiny ( ccum ) ) then
    ccum = 0.0D+00
  end if

  return
end subroutine cumnor      

subroutine r8_swap ( x, y )

!*****************************************************************************80
!
!! R8_SWAP swaps two R8 values.
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none

  real(8) :: x
  real(8) :: y
  real(8) :: z

  z = x
  x = y
  y = z

  return
end subroutine r8_swap

END MODULE STATISTICS
